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PREFACE 


This  Memorandum  represents  a  minor  effort  under  the  Project 
RAND  program  for  general  geophysics.  It  has  been  noted  that  data 
from  the  ROBINs  (rocket  balloon  instruments)  appear  inconsistent 
and  that  these  falling  spheres  appear  to  collapse  more  often  than 
would  seem  reasonable.  This  Memorandum  may  offer  a  significant 
insight  into  why  some  of  the  data  from  the  ROBIN  experiments  are  lost 
and  how  some  may  be  misconstrued.  It  should  be  of  general  interest 
to  groups  studying  the  characteristics  of  the  mesosphere  or  operating 
rockets  in  and  through  that  region. 
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ABSTRACT 

In  order  to  study  possible  perturbations  of  the  ROBIN  density 
calculations,  a  consistent  set  of  velocities  was  computed  between 
70  km  and  38  km  for  a  standard  atmosphere  using  simplified  fall-velocity 
equations.  Consistency  was  achieved  by  iterating  the  calculations 
until  variations  in  the  drag  coefficient  produced  insignificant  changes 
in  the  fall  velocities.  Fall  velocities  were  then  computed  using  an 
equation  that  included  an  arbitrary  distribution  of  vertical  atmospheric 
motions.  The  resulting  fall  velocities  were  used  to  compute  air  density, 
assuming,  as  is  done  in  the  ROBIN  calculations,  that  there  is  no 
vertical  air  motion.  These  computations  are  used  to  demonstrate  the 
error  in  density  resulting  from  the  neglect  of  the  vertical-motion  term. 
The  effect  of  the  vertical  motion  on  the  iterative  process  used  to 
obtain  the  drag  coefficients  and  on  the  process  used  to  check  the 
balloon's  sphericity  is  also  discussed. 
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I .  INTRODUCTION 


The  ROBIN  (Rocket  Balloon  Instrument)  is  a  spheroid  of  1/4-mil 
Mylar,  1  meter  in  diameter,  having  a  metalized  Mylar  corner  plate  used 
as  a  reflector,  which  when  ejected  from  a  rocket  at  some  altitude, 
inflates  and  is  tracked  by  radar  as  it  falls.  With  the  information 
recorded  by  the  radar,  the  ROBIN  becomes  a  tool  for  determining 
atmospheric  density.  Because  of  seemingly  inconsistent  density  data 
from  ROBIN  experiments  (Ref,  1,  Fig.  5),  this  study  was  undertaken  to 
examine  vertical  winds  as  a  possible  cause  of  these  inconsistencies 
in  the  computed  densities. 

One  of  the  first  problems  to  be  solved  is  that  of  determining 
proper  drag  coefficients.  Since  drag  coefficients  are  not  known  as  a 
function  of  altitude  and,  in  reality,  cannot  be  determined  without 
knowledge  of  the  Mach  and  Reynolds  numbers,  it  is  necessary  to  assume 
their  values  and  observe  what  errors  may  ensue  by  back  calculation.  As 
a  result  of  his  work  with  the  ROBIN,  Engler  (Ref.  2,  p.  93)  presents 
a  table  of  drag  coefficients  as  a  function  of  Mach  and  Reynolds  numbers, 
which  provides  the  values  necessary  for  the  initial  assumptions.  If 
the  calculations  are  iterated,  variations  in  drag  coefficients  will  be 
reduced  to  the  point  of  producing  insignificant  changes  in  the  velocity. 

The  introduction  of  vertical  winds  into  the  problem  produces  interesting 
changes  in  velocity  measurements  and  density  calculations.  It  is  not 
possible,  from  ROBIN  measurements,  to  distinguish  between  changes  in 
density  and  vertical  wind  perturbations;  we  will  discuss  this  dependency 
of  calculated-density  variations  upon  vertical-wind  magnitudes.  Also, 
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II.  DEVELOPMENT  OF  A  CONSISTENT  SET  OF  P  VS,  w  VALUES 


Using  the  coordinate  system  in  which  altitude,  z,  is  positive 

(3) 

downward,  the  equations  of  motion  of  Lally  and  LevitonN  7  become 


dw 


dt 


in  which 


g  -  aT 


2  2 
Cdttt  pw 


2m 


where  g  is  gravity,  CD  is  the  drag  coefficient,  r  is  the  radius,  m  is 

the  mass  of  the  balloon,  p  is  the  density,  and  w  denotes  the  balloon's 
* 

vertical  velocity.  Thus 


dw  KC_pw 

—  -  g - 

dt  2 


a) 


where  K  *  to:  /m  is  a  constant.  Since 


dw  dw  d(w  /2) 
dt  dz  dz 


vertical  acceleration  can  be  expressed  as  a  function  of  altitude  as 


d(w2/2) 


dz 


g  - 


KCDpv 


(2) 


★In  this  case,  the  air  is  assumed  to  be  still,  and  the  balloon's 
velocity  relative  to  the  earth  is  the  same  as  the  balloon's  velocity 
relative  to  the  air. 
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any  variation  in  drag  coefficient  will  result  in  a  variation  in 
computed  density;  therefore,  as  will  be  shown,  vertical  winds,  drag 
coefficients,  and  densities  are  quite  closely  related. 

Using  the  simple  equation  defining  the  fall  rate  of  a  ROBIN  as 
determined  by  Lally  and  Levitonv  7  and  the  equations  defining  a 
standard  atmosphere  (Ref.  4,  pp.  2,6),  we  have  incorporated  vertical 
winds,  as  proposed  by  Woodbridge^  in  the  calculation  of  the  balloon's 
vertical  velocity.  The  hypothesis  of  the  presence  of  vertical  winds  will 
be  further  illustrated  by  the  use  of  data  from  ROBIN  flights  (Ref.  2, 
pp.  3-4)  to  calculate  these  vertical  perturbations. 

A  discussion  of  the  \ -check  as  proposed  by  Engler  (Ref.  2,  pp.  11-14) 
is  included  to  point  out  the  difficulty  of  differentiating  between 
fall  rates  resulting  from  a  vertical  wind  and  those  resulting  from 
malfunction  of  the  balloon. 

The  IBM- 7090  was  used  to  perform  the  numerical  calculations  to 
determine  the  balloon's  velocity.  An  appendix  is  included  to  illustrate 
the  mechanisms  of  the  FORTRAN  program. 


I* 
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2 

By  defining  q  =  w  /2,  Eq.  (2)  becomes 

dq 

—  =  g  "  KCDpq  ,  (3N 

dz 

which,  by  numerical  integration  (see  Appendix),  yields  the  vertical 
velocity  of  the  balloon. 

Standard  density  and  temperature  are  calculated  by  use  of  a  con¬ 
sistent  set  of  equations,  i.e.,  those  applying  to  the  standard  atmosphere, ^ 
which  are  (in  a  coordinate  system  with  z  positive  downward)  as  follows: 

dP 

—  *  Pg 
dz 

P  =  PRT  (4) 

T  =  Tb  -  v(z  -  *b) 

with  mks  units  where 

-2 

g  =  gravity  (m  sec  ) , 

-3 

p  =  density  (kg  m  ), 

-1  -9 

P  -  atmospheric  pressure  (kg  m  sec  ), 

2  o—l—l 

R  -  universal  gas  constant,  2.87  x  l(r  joules  (°K)  kg  , 

T  -  temperature  (°K) , 

T,  =  value  of  T  at  altitude  z.  , 
d  b 

y  =  temperature  gradient  (°K  m"*^). 

Using  Eqs ,  (4)  to  solve  for  density,  p,  we  find 

dP 

—  =  Pg  , 
dz 


i 
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dP  dp  dT 

—  *  RT  —  +  Rp  —  , 

dz  dz  dz 

dT 

—  *  -Y  . 

dz 


So, 


and 


dp 

RT - RpY  -  Pg  -  0, 

dz 


1  dp  g  +  Ry 


p  dz  RT 


or 


dinp  g  +  Ry 

dz  R[Tb  -  y(z  "  zb)l 


Upon  integration  of  Eq.  (5)  from  z^  to  z,  where  z^  is  the 
the  base  of  the  interval, 


P 


(1  +  g/Ry) 


“  Y (z"zb) 


for  y  t  0. 

For  y  “  0,  Eq.  (5)  becomes 


dXnp  g 

d2  RTb 


yielding 


pb  exp 


8(z 


RT, 


c 


(5) 

altitude  at 


(6) 


(7) 
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The  equation  used  to  calculate  temperature  is  simply 

T  -  *b  -  V(«  -  zb> 

as  noted  above. 

Drag  coefficients  present  a  problem  in  these  calculations  for 
velocity  since  they  are  not  expressible  as  a  function  of  altitude. 

Because  the  drag  coefficients  are  actually  unknown  until  a  velocity  is 
determined  and  Mach  and  Reynolds  numbers  can  be  found,  it  is  necessary 
to  assume  a  drag  coefficient,  C^,  to  compute  a  velocity,  then  by 
computing  the  corresponding  Mach  and  Reynolds  numbers,  determine  a 
new  drag  coefficient,  C^,  to  be  compared  with  the  original  assumption 
as  a  check. 

We  made  the  original,  arbitrary  assumption  of  a  set  of  values, 
using  the  drag  coefficients  of  Engler  (Ref.  2,  p.  93)  as  a  basis;  these 
are  shown  in  Fig,  1  as  C^.  With  these,  we  determined  corresponding 
velocities  and  Mach  and  Reynolds  numbers,  as  discussed  above,  to  obtain 
a  more  refined  drag  coefficient  for  the  problem;  the  results  of  this 
process  to  obtain  a  new  drag  coefficient,  C^,  are  also  shown  in  Fig.  1. 

Then  is  the  drag-coefficient  value  used  in  all  subsequent  calculations. 

Upon  comparison  of  the  original  and  the  refined  drag  coefficients,  it 
can  be  seen  that  varies  from  by  07.  to  2.5%.  Thus,  by  refining 
the  drag  coefficient  values  only  once,  a  more  accurate  approximation 
of  velocity  can  be  achieved. 

Let  us  consider  the  results  of  the  refined  drag  coefficients  on  the 
subsequent  velocity  calculations.  A  range  of  deviation  of  drag- coefficient 
values  of  0%  to  2.5%  results  in  a  range  of  deviation  for  velocity 
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Fig.1  —  Drag  coefficient  originally  assumed  and 
refined  drag  coefficient 
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calculations  of  0%  to  1.3%.  So  it  can  be  seen  that  absolute  accuracy 
in  drag-coefficient  values  is  not  required  for  approximations  of  velocity 
within  a  minimal  range  of  error  in  measurement. 
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III.  INTRODUCTION  OF  VERTICAL  WIND 

Using  the  vertical-wind  data  plotted  by  Woodbridge^  as  a 
,  * 

basis,  we  constructed  a  smooth  arbitrary  curve  over  the  30-km  interval 
between  40  km  and  70  km  in  which  the  initial  and  terminal  vertical 
wind  velocity,  defined  as  u,  was  zero.  From  this  curve,  values  of  u, 
selected  at  given  intervals  as  a  function  of  altitude,  were  entered 
into  the  calculations  for  the  vertical  velocity  of  the  balloon,  w. 

Upon  consideration  of  a  vertical  wind  not  equal  to  zero,  the 
equation  of  motion  (1)  becomes,  after  replacing  w  with  (w-u)  in  the 
drag  term, 

dw 

—  -  g  -  KC'p[(w-u)z/2  J  .  (8) 

dt  U 

From  Eq.  (3),  we  find 
dq 

-  =  g  -  KC -p 
dz 

where  z  is  positive  downward;  the  drag  force  is  now  proportional  to  the 
relative  vertical  wind,  and  the  balloon* s  vertical  acceleration  is 
relative  to  the  earth.  Since,  over  each  subinterval  of  integration, 
g,  K,  C^,  p,  and  u  are  given  as  constant,  this  equation  is  easily 
solved  for  q,  and  thus,  w  is  determined. 

The  vertical  velocity  of  the  balloon  as  determined  by  Eq.  (9) 
is  now  a  function  of  a  vertical  wind,  u,  as  well  as  being  a  function 

^Although  the  values  from  Woodbridge  are  used  as  a  basis  for  defining 
a  possible  vertical  wind  velocity  distribution,  we  do  not  mean  to  imply 
that  his  values  are  accepted  as  correct. 
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of  altitude,  z .  We  shall  designate  this  velocity  as  and  show  that 
when  it  is  used  to  calculate  density  according  to  Eq.  (12),  below,  in 
an  atmosphere  containing  a  vertical  wind,  an  erroneous  density,  p^, 
results . 

By  assuming  that  altitude  is  positive  downward  and  taking  the 
differential  equation 

d(v2/2)  , 

-  =  g  -  KCpP  (*»  /  2)  ,  (10) 

dz 


by  defining  the  balloon’s  motion  with  no  vertical  wind,  by  making  the 
2 

substitution  for  w  /2,  and  by  solving  for  p  by  taking  finite  differences, 


we  find  for  Eq,  (10) 

1  Aq 

p  =  -  g  +  — 

KCi  q  L  Az. 


(ii) 


When  u  f  0,  a  similar  form  can  be  found: 


P 


u. 

l 


“d^i  L 


gi  + 


(qi-l  -  qi+l 
Az 


(12) 


where  i  denotes  the  sub  interval  of  evaluation.  Using  this  method  of 
finite  differences,  the  balloon’s  vertical  velocity,  including  a 
vertical  wind,  determines  a  density,  p^.  It  should  be  mentioned 
in  passing  that,  even  though  errors  are  small,  there  is  some  error  in 
the  results  when  the  finite-difference  method  of  back  calculation  is 
used,  as  compared  with  the  results  supplied  by  the  IBM-7090  computer. 

The  results  of  Eqs.  (9)  and  (12),  i.e.,  the  results  of  the 
calculations  for  velocity  and  density  assuming  nonzero  vertical  air 
motion,  can  be  seen  in  Fig.  2.  Figure  2a  illustrates  the  arbitrarily 
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The  effect  of  ossumed  vertical  wind  velocities  on  velocity  of  falling 
ROBIN,  and  the  resulting  miscalculation  of  density 
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chosen  variation  in  the  vertical  wind;  Fig.  2b  shows  the  effect  of 
this  vertical  air  motion  on  the  vertical  motion  of  the  balloon;  Fig.  2c 
gives  the  erroneous  density  corresponding  to  the  shown  velocity 
variation  which,  in  turn,  is  a  result  of  the  given  vertical  wind.  As 
should  be  expected,  an  upward  air  thrust  results  in  a  lesser  balloon 
velocity  and  a  greater  calculated  density;  downward  air  motion 
results  in  the  opposite  —  greater  balloon  velocity  and  a  smaller 
calculated  density. 
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IV.  CALCULATION  OF  PERCENTAGE  ERROR  DUE  TO  VERTICAL  WINDS 


Since,  when  a  nonzero  vertical  wind  exists,  an  erroneous  calculated 
density  directly  results  as  discussed  in  Section  III,  it  is  necessary  to 
be  able  to  make  an  approximation  as  to  what  magnitude  relative  error 
should  be  expected.  Here  we  will  develop  an  expression  for  the 
approximation  of  relative  errors  in  density  calculations  as  a  function 
of  a  vertical  air  motion,  u,  and  thus,  show  the  dependence  of  density 
calculations  on  this  vertical  motion.  Equation  (1),  u  =  0,  becomes 

d(w2/2)  Spstdw2 

-  =  g - , 

dz  2 

and  Eq.  (8),  u  f  0,  becomes 

d(w2/ 2)  K(^pu(w-u)2 

-  "  g - 

dz  2 


Therefore, 


or 


2(^  -  8>/<-K Cj)  Pstd“2 

'  8>/<-%>  pu("‘u>2 


Pu(w-u)2 


and 


^std 

pu 


(w-u)2 


(13) 
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By  expanding  and  simplifying  Eq.  (13), 


and 


p  w  w 


2u 

w 


(u/w)2  . 


(14) 


Therefore,  the  relative  error  in  the  determined  density  may  be  expressed 

2  2 

as  Eq.  (14).  Since  u  «  w,  the  term  (u  /w  )  in  Eq.  (14)  was  omitted; 
its  contribution  to  the  determination  of  relative  error  in  density 
calculation  is  negligible.  The  approximation  of  the  relative  error 
in  the  calculated  density  is  made  primarily  by  (2u/w) .  Thus 


Pu  " 


std 


P 


u 


2u 

w 


Let  us  now  assume  that  we  should  expect  a  density-calculation 
error  of  2%  to  10%;  we  wish  to  know  the  vertical  wind  velocity  that 
would  cause  a  relative  error  in  this  range.  Since  the  relative 
error  in  density  calculation  is  approximated  by  2u/w,  we  can  say, 
for  2%  error  in  density,  that  0.02  tv  2u/w,  or  u  tv  O.Olw;  for  a  10% 
error  in  density,  u  ^  0.05w.  Figure  3  illustrates  the  balloon's  vertical 
velocity  with  no  vertical  wind  perturbation  and  the  vertical  winds 
that  will  cause  2%  and  10%  density-calculation  errors.  It  can  be 
seen  here  that,  in  the  hypothesized  model,  density-determination 
accuracy  and  vertical-wind  magnitude  are  inseparable. 
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.  THE  EFFECT  OF  VERTICAL  WINDS  ON  THE  ITERATION  FORMULA. 


Vertical  winds  will  introduce  errors  into  the  determination  of 
drag  coefficients  which  will,  in  turn,  cause  an  error  in  calculated 
density  additional  to  that  caused  by  vertical  wind  directly.  The 
procedure  below  was  followed  to  show  these  variations  in  drag-coefficient 
and  density  calculations.  With  the  values  for  C^  and  p  and  the  resulting 
computer  calculations  for  velocity,  w^,  Mach  and  Reynolds  numbers 
were  determined  from 

w  2pw  r(T4-S) 

m  =  -2-  p  -  u  — 

KT^  ’  e'  B3'2 

where  K  =  20*05  m  sec  \°K)  2,  $  =  1.458  x  10  ^  kg  m  ^  sec  ^(°K)  \ 

S  =  110. 4°K,  T  is  the  atmospheric  temperature,  and  r  is  the  radius 
of  the  balloon;  M  and  Re  then  determine  a  drag  coefficient  value,  Cjj, 
by  use  of  Engler's  table  (Ref,  2,  p.  93).  With  CjJ,  a  new  density 
value,  p*  is  found  from 


which  is  similar  to  Eq.  (11).  As  can  be  seen,  the  variation  in  the 
density  is  inversely  related  to  the  variation  in  the  drag  coefficient. 
In  Table  1,  a  few  examples  of  this  are  shown.  Whenever  indeterminable 
vertical  winds  exist  there  will  be  errors  in  drag-coefficient  deter¬ 
minations  from  the  ROBIN  data  and,  therefore,  errors  in  the  relative 
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density,  Thus  it  is  evident  that  in  the  presence  of  vertical  winds, 
density  measurements  will  not  necessarily  converge  to  a  correct  value. 

Table  1 

VARIATION  OF  DRAG -COEFFICIENT  AND  DENSITY  CALCULATIONS 


P  (kg  m  3) 


p'(kg  m'3) 


km 

.473  .468 

1.657  x  lo"4  1.824  x  lo"4 

66 

.455  .446 

2.717  2.771 

62 

.435  .426 

4.325  4.416 

58 

.426  .427 

6.707  x  lo’4  6.694  x  lo"4 

54 

.430  .434 

1.045  x  lo"3  1.036  x  10'3 

50 

.437  .443 

1.671  x  lo"3  1.648  x  10-3 

46 
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VI.  THE  EFFECT  OF  VERTICAL  WINDS  ON  THE  X  -CHECK 

An  approximation  of  the  X -check  is  used  by  Engler  (Ref.  2, 
pp.  11-14)  as  a  method  for  determining  whether  or  not  the  balloon  is 
spherical.  With  a  properly  inflated  balloon  and  no  vertical  wind,  the 
value  of  a  calculated  parameter,  X,  can  be  shown  to  be  limited  to  a 
specific  range.  The  existence  of  a  vertical  wind  acting  on  a  ROBIN 
might  produce  the  same  effect  during  a  X -check  as  a  balloon  malfunction. 

The  value  X  in  Engler' s  method  is  approximated  by  Xq,  which  is 
defined  by 

2(dw/dt)  2z 

Xrt  -  - r  -  - r 

(dz/dt)Z  (i)2 

At  the  maxima  and  minima  of  dw/dt,  XQ  is  computed  by  Engler,  and 

if  X  >  2  x  10  ^  m  *,  or  X  <  0.5  x  10  ^  n  \  it  is  assumed  that  the 

o  o 

balloon  is  not  rigid,  and  density  calculations  stop. 

By  definition, 

dw  dw  d(w^/2) 

z  ■  —  *  w  —  *  ■  ■ 

dt  dz  dz 

Therefore,  since  z  ■  w, 

S  -  g  -  KpCD  [(£)2/2  J  , 

and 

(i)2/ 2  “  q  ; 

80*  z  g  -  KpCDq 


q 


q 
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or 

Xo  .  g/,  -  k»cd 

Table  2  below  shows  the  results  of  the  X -check  on  the  maxima  and 
minima  of  the  acceleration  in  our  model  ROBIN  experiment  as  determined 
by  the  IBM- 7090. 


Table  2 

RESULTS  OF  THE  X -CHECK  ON  THE  ACCELERATION  OF  A 
ROBIN  RESPONDING  TO  A  NONZERO  VERTICAL  AIR  MOTION 


Accel¬ 

eration 

u 

Cm/ sec) 

w 

u 

(b/ sec) 

g 

2 

(m/sec  ) 

cp 

P 

_£kgZmi)  __ 

q 

(m/aec)^ 

fin”1! 

Max 

-5. 

117.14 

9.63 

.426 

5 .402xlo"4 

6.861xl03 

.23x10'^ 

Min 

+4. 

83.69 

9.64 

.426 

5 . 705xl(f 4 

7 .004xl03 

-.82x10”* 

As  can  be  seen,  the  X -check  at  both  the  maximum  and  minimun  accelerations 
yields  a  ^  <  0.5  x  10  \  and  density  measurements  should,  according 

to  the  criteria,  cease.  Lambda  measurements  at  the  maximum  and  minimum 
accelerations  would  normally  be  considered  to  be  valid  since,  at  this 
point,  the  ROBIN  has  passed  its  peak  velocity.  Here  we  can  see  an 
example  of  a  vertical  wind  being  mistaken  for  a  malfunction  in  the 


balloon 
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VII.  CONCLUSIONS 


The  process  of  deducing  the  atmospheric  density  from  the  actual 
tracking  of  ROBIN  spheres  is,  of  necessity,  based  on  the  assumption  of 
zero  vertical  wind.  By  constructing  a  consistent  set  of  hypothetical 
vertical  sphere  velocities  based  on  a  standard  atmosphere,  and  a  given 
set  of  vertical  winds,  we  have  been  able  to  illustrate  the  effect  of 
vertical  winds  on  the  computed  densities.  Our  results  indicate 

(1)  that  vertical  winds  of  reasonable  velocity  can  directly 
increase  the  errors  of  ROBIN-derived  densities  beyond  the 
published  error; 

(2)  that  these  vertical  winds  affect  the  iteration  technique 
for  determining  drag  coefficients  to  a  degree  that  casts 
doubt  on  the  validity  of  the  drag  estimates; 

(3)  that  they  might  also  cause  accelerations  that  could  be 
mistaken  for  malfunctions  of  the  balloons. 

These  indications  suggest  that  the  falling-sphere  technique  needs 
a  critical  re-examination.  Perhaps  a  relaxation  of  the  stated 
accuracy  and  stated  resolution  of  the  measurements  of  the  spheres  would 
produce  more  useful  and  more  nearly  consistent  results. 


Appendix 
COMPUTER  PROGRAM 


Hie  FORTRAN  program  discussed  here  was  written  for  the  IBM- 7090 
to  provide  a  step-wise  solution  to  the  first-order,  differential 
equation,  Eq .  (1),  defining  the  fall  rate  of  the  ROBIN.  Although  a 
block  diagram  is  included  (Fig.  4),  some  discussion  is  necessary. 

The  solutions  were  obtained  by  means  of  a  Runge-Kutta  fourth- 
order  method,  modified  by  Blum.^  Conte  and  Titus ^  discuss  the 
integration  technique  that  is  incorporated  in  RAND  routine  (FORTRAN) 
RS-X013;  it  is  also  available  in  the  SQUOZE  deck  as  X004.  The  calling 
sequence  for  the  set-up  of  the  integration  routine,  which  must  be 
used  before  the  integration  is  started  and  once  per  solution,  is 
CALL  INT  (T,  N,  Kl,  A2,  A3,  A4,  A5,  A6,  A7) 
where  T,  a  region  of  at  least  12N^3,  is  defined  as  follows: 

T(l) :  set  up  by  CALL  INT 

T(2):  value  of  independent  variable,  z 

T(3):  value  of  step-size,  Az 

T(4):  value  of  dependent  variable,  q,  supplied  initially  by  user 
T(5):  value  of  derivative,  dq/dz,  stored  by  DAUX. 

The  calling  sequence  for  integrating  one  step  is  CALL  INTM  and  requires 
no  arguments.  Over  each  increment  of  integration,  Az,  INTM  calls  up 
DAUX  four  times  to  evaluate  the  derivative  dq/dz.  The  parameters  of 
the  integration  routine  are  the  following: 

N  :  number  of  equations 

Kl  :  option  word  for  determination  of  integration  mode  (0  applies 

to  the  Adams -Moulton  variable  step-size  mode,  1  applies  to  the 
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MAIN  ROUTINE 


SUWOUTI  NE  OAUX 


Fig.  4 —  Block  diogrom  of  FORTRAN  program  to  study 

ROBIN  fall  speeds 
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Runge-Kutta-Blum  mode  used  for  this  problem,  and  2  applies  to 
the  Adams-Moulton  fixed  step-size  mode.) 

A2,  A3,  AA,  A5,  A6,  A7  : 

apply  only  to  the  Adams-Moulton  variable  step-size  mode. 

A  sample  (Fig.  5)  of  the  input  formats,  along  with  a  table  that 
relates  FORTRAN  nomenclature  to  the  text,  is  included  to  illustrate 
the  manner  in  which  the  data  are  read,  as  listed  in  the  main  routine 
and  subroutine  (Fig.  6).  A  sample  of  the  output  format  is  also  included 
(Fig.  7). 

A  comment  should  be  made  regarding  the  technique  used  in  the 
FORTRAN  program  to  determine  the  fall  speeds  of  the  ROBIN.  Each  given 
2-km  interval  is  divided  into  twenty  100-meter  subintervals.  Values 
for  C^,  drag  coefficient,  g,  gravity,  and  u,  vertical  wind,  are  read 
as  data  and  held  constant  over  each  given  2-km  interval  (controlled 
by  I);  p,  density,  and  T^,  temperature  at  the  base  of  an  interval, 
having  only  initial  values  supplied,  are  calculated  over  every 
100-meter  subinterval  (controlled  by  K) .  Values  of  p  and  Tfe  are 
assumed  to  change  very  little  over  100  meters  and  therefore,  are 
calculated  outside  of  DAUX  and  held  constant  over  each  subinterval, 
or  increment,  of  evaluation. 
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i 


Somple  imput  format 
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Table  3 

DEFINITION  OF  FORTRAN  NOMENCLATURE 


FORTRAN  Program 


R 

RHO 

RSTD 

TB 

QZ 

WIND 

WZ 


r,  radius  of  the  balloon 

p,  density 

R,  universal  constant 

T^,  temperature  at  base  of  given  interval 
2 

q,  defined  as  w  /2 
u,  vertical  wind 

w,  velocity  of  the  balloon 


Z 


altitude 
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COMMON  T.N  .  I  *  K  *  P  .CONS  T  *RHO  *  Z  •  COX  *  GX  .GAMMA X  .  ToX . VW I NDX . DZ • Q 
DIMENSION  T<5  .0  ) .RHO( 350) t CDU 6 >  tG O A ) . VW I  NO (1 6 I . TB I  1 6 ) .GAMMA ( 16 ) 

1 tCZ(35G) .WZl 16) .Z( 16) 

INPUT  CODE  NUMBER  FUR  FORTRAN  INTEGRATION  ROUTINE  (Kilt  AND  THE  NUMBER 
OF  INTERVALS  OF  INTEGRATION  (IT). 

READ  1  *  K1 *  IT 

]  F  Ci  Li  7  k  2  i  1  2  > 

PRINT  CODE  NJMotk.  NUkBOk  Of  INTERVAL  INI  t ORATION.  AND  INTEGRATION 
MODE  USED. 

PRINT  10.  Kit  IT 
10  FORMA  T ( 2 1  1 0 ) 

PRINT  20 

20  i-OKMAT  ( 24H0RUNGE-KUTTA  INTEGRATION) 

2  FORMA T ( 6E 1 2 • 8  ) 

3  FORMA T ( 3E12.8 ) 

INPUT  AND  PRINT  INTEGRATION  ROUTINE  FORMAT • 

READ  2*  A2*  A3.  A4*  A5  .  A6 *  A7 
PRINT  30.  A2*A3  .A4.AS.A6.A7 

30  FORMAT (1CCH  EbAK  M  A 

1  nZMAX  DZMIN  BETA/ 1H  6E18.6) 

PRINT  40 

40  FORMAT  <8H0DZ=2  KM) 

INPUT  DRAG  COEFFICIENT.  GRAVITY.  AND  VERTICAL  WIND  WHICH  WILL  dE  CON¬ 
STANT  OVER  EACH  INTERVAL  BUT  WILL  CHANGE  EVERY  2KM. 

ALTITUDE  IS  POSITIVE  DOWNWARD. 

RFAD  3.  ( CD( J ) tG( J) tVWINM J ) •  J=1*IT) 

INPUT  THE  MASS  AND  RADIUS  OF  THE  ROBIN.  CALCULATE  CONSTANT.  K.  PRINT 
HEADING  AND  OIJTPJT. 

READ  3.  R.  oKASS 
PI *3. 1 415? 

CONST =PI*(R**2 ) /BMASS 
PRINT  SO.  R.DKASS. CONST 

50  FORMAT (65HO  3  BM ASS  CONST 

1  /3E20.3I 

SET  UP  INITIAL  CONSTANTS  AND  CONDITIONS  FOR  INTEGRATION. 

QZ(1 ) =C. 

DZ  =1^U 

RH0( 1  I =9.906/ ( 1  0***5) 

TB(1) =216.36 
T(2)=-. 

T~(~3  >=DZ 
T  ( 4 ) =OZ ( 1 ) 

T(5)=9,59 

CALL  SET-UP  OF  INTEGRATION  RC JT I NT 

CALL  INT ( T  *1 *K1 .A2.A3.A4.A5  .A6.A7) 

dEG*N  DO  LOOP  CONTROLLING  NLT'BlK  OF  I  NT f° VALE  OF  INTEGRATION. 

DO  666  I  =  1  .  I  T 
1  =  1 

SET  VALUES  OF  DRAG.  GRAVITY.  AND  VERTICAL  WIND  FOR  2 K V.  INTERVAL. 

CDX  =  C  D ( I  ) 

GX  =  G(  I  ) 

V  W  I  N  D  X  =  V  a  I  N  D  (  I  ) 

test  altitude  to  determine  valul  of  tlmpepature  gradient,  gammax. 

IF(T(2)-23O50.)  51.667*32 

51  IFC17C50.-TC2) )  53.667.32 

52  GAMMAX= - .0039 

GO  TO  55 

53  GAMMA X  =  j. 

,v:  TO  05 

54  G A  Mf/ AX  =  •  00  30 

53  CONTINUE 

DEFINE  Q.  RSTD (UNIVERSAL  CONSTANT).  P  TC  BE  USED  LATER. 

Q=AB5F‘  ( GAMMAX  ) 

RSTD  =  28  7.0 

H=1 •  +  ( GX / ( R  o  T  D* GAMMAX  )  ) 


Fig.  6  —  Program  listing 


I 
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print  HEADING. 

PRINT  60 

6C  FORMAT ( 60H1  RHO  T3  GZ 

1  ) 

BEGIN  DO  LOOP  WHICH  CONTROLS  CALCULATION  CF  DENSITY.  TEMPERATURE.  AND 
1  (  W  ft  ft  2  )  /  2  )  OVER  2K.M  INTERVAL. 

DO  66  <=1.20 

K  =  < 

INTEGRATE  ONE  ALTITUDE  INCREMENT. 

CALL  INTM 

PRINT  70*  KHCiKl .TB(<) *T(4) 

7  U  FORMAT ( 1H03E20.B) 

CALCULATE  TB  FOR  NEXT  INTERVAL. 

TbU+l  )  =Tb(  <)-(GAMMAX*DZ) 

TEST  GAMMAX  AND  CALCULATE  DENSITY  ACCORDINGLY. 

1MC)  667.61.6? 

61  Rn0(K+l ) =RhO(< ) *< EXFF ( < GY *01 ) / ( RSTD*T3 ( K ) m 

GO  TO  66 

6  2  RHO ( <  + 1 ) =  RhO ( < )* ( TB( <) / ( TQtK+1 )  >  )**P 

65  CONTINUE 

66  CONTINUE 

3E I  CALCULATED  VALUES  of  density  and  temperature. 

RHOX=RHO( 20) 

TbX*Tu ( 20 ) 

TEST  NJMdLk  OF  INTEORaTlD  I  NT  ckVALd.  IF  INCOMPLETE*  CONTINUE.  IF  COM¬ 
PLETE.  PRINT  FINAL  ANSWERS  AND  STOP. 

I F ( I  -  I T )  67.75  *75 

67  RHQ( 1 ) =RHC< 21  ) 

TB<  1  J*TB  t  2U 

75  CONTINUE 

WZ( t ) =SQRTf  (2.*T14)  ) 

Z( I ) =  7 . * ( lu.**4)-T<2  ) 

PRINT  60 

80  FORMAT!/////) 

PRINT  9 x. 

90  FORMAT llHP tlCX»lHZ . 1 6 K ,2 hCD *  1 6 X , 6HGAKMA *  1 7X  * 1HG  *  1 7X  1 5HVW 1  ND/ 

1  1H  *9X  .  2hOZ  >  l  ttX  .2Ha‘Z*  18X  .  2HTb  *  1  8X  .5HRH0  > 

PRINT  10  ,.Z(I)  *CDX  *  GAMMA  X  *GX  .  V.*.  I  NDX  .  T  (  4  )  *  v  L  (  I  )  .TbX  ,RHOX 
10 u  F0RMAT(1H05E2C.8/(1H  4E20.8)) 

666  CONTINUE 
GO  TO  1000 

667  PRINT  668 

668  FORMAT  U9H1S0ME  THING  IS  WRONG) 

100 u  CONTINUE 

CALL  EXIT 
END 


SUBROUTINE  DAUX 

COMMON  T  .N*  I  *<.P. CONST  .RHO  *  l  ♦  COX  •  GX  *GA  ”MA  X  .  T3X  *  V.V  I  NDX  .  DZ  *  0 
DIMENSION  T(5vO) .RHO (350) *CD ( 1 6 ) . G ( 1 6 ) . VW I ND ( 1 6 ) » T 9 ( 1 6 ) .GAMMA ( 16) 
1  .QZ  (  350  .  WZ  ( 16  )  *Z  ( 16  ) 

ROo IN’S  EQUATION  OF  MOTION  INCLUDING  VERTICAL  WIND  TERM. 

T ( 5 ) =GX-(CONST*RHGl<)*CDX)*( T ( 4 ) -VW I NDX* ( SO^TFl 2 .*T<  4 )  )  ) 

1  +( ( VWlNDX)**2)/2. ) 

RETURN 

END 


Fig.  6  ( cont) — Program  listing 
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EXECUTION 

1  16 

RUNGE-KUTTA  INTEGRATION 
ERAR  H  A 

0.  0.  0 

DZMAX 

.  0. 

DZMIN 

0* 

BET 

0. 

DZ-2  KF 

R 

0. 5CC00CCCE  00 

BFASS 

0*  13000000E-00 

CCNST 

0.60415192E  01 

RHC 

TB 

cz 

0.59059<5S8F-0<> 

0 • 216  36000  E 

03 

0  •  78807286E 

03 

0 • IC04  193  56-03 

0 • 2 167  5000E 

03 

0*  17028362E 

04 

0. 10179987F-03 

0*21  7 14000E 

03 

0*258474886 

04 

0,  103186  71 E-0  3 

0.21752999E 

0*3 

C* 3434111 3E 

04 

0. 10459504E-03 

0,21 79 1999E 

03 

0.42512459E 

04 

0,  10602C02  E-0  3 

0 • 2 18  30999E 

03 

0  •  503  64958E 

04 

0,1 0746 1 82  E-03 

0*21 969999E 

03 

0 . 5  75  02249E 

04 

0, 10892062E-03 

0.21908999E 

03 

0.651  281 66E 

04 

0,1 1039656E-0  3 

0.21947999E 

03 

0 • 72C4674 1 E 

04 

0, 1 1168982E-03 

0.21986999E 

03 

0 • 786  62 1 92E 

04 

0,11 340059F-03 

0 • 2202  5999E 

03 

0*  84978916E 

04 

0,1 14929C  3E— C  3 

0  *22064  9996 

03 

0  •  9 1 C  0 1490E 

04 

0*1164  7531  E-0  3 

0 • 22 103999E 

03 

0*967 34650E 

04 

0*1 180  3962E-03 

0  •  22 142999E 

03 

0*102 18  330E 

05 

0*1 19622 1 3E-03 

0 *  22 18  1999E 

03 

0*  107  35248E 

05 

0*1212230 1 E-03 

0 • 2222C998E 

03 

0. 11224740E 

05 

0* 12284246 E-Cj 

0 • 22259990  E 

03 

0.  U637339E 

05 

0* 12448 06 5E -03 

0 • 22298998E 

03 

0*  121 2  3590E 

05 

0*1261 3  776E-03 

0 • 22  337998  E 

03 

0  •  125  34051 E 

05 

0 • 12781 399F— 03 

0*2?  376999E 

03 

0. 12919292E 

05 

l  CO 

CZ  W  Z 

GAFFAt 

Tti 

G 

RHO 

VWIND 

0.68CC0C00E  05  0.5849S999E  00  -0. 38S99999E- 

0.12919292E  05  0.  1607A385E  03  0.22376999E 

*02  0 • 95899999E  01  0. 

03  0*12781 399E-03 

Fig.  7 — Sample  output  format 
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